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The electronic properties of graphene are influenced by both geometric confinement and strain. 

We study the electronic structure of in-plane bent graphene nanoribbons, systems where confinement 
and strain are combined. To understand its electronic properties, we develop a tight-binding model 
that has a small computational cost and is based on exponentially decaying hopping and overlap 
parameters. Using this model, we show that the edge states in zigzag graphene nanoribbons are 
sensitive to bending and develop an effective dispersion that can be described by a one-dimensional 
atomic chain model. Because the velocity of the electrons at the edge is proportional to the slope 
of the dispersion, the edge states become gradually delocalized upon increasing the strength of 
bending. 


I. INTRODUCTION 

Many of graphene’s remarkable features stem from two 
facts. The first is that its low energy quasiparticles are 
linearly dispersive and can be effectively described as 
Dirac fermions® the second is that graphene is a two- 
dimensional ultrathin membrane that holds promises to 
revolutionize the current nanotechnology.® In addition, 
this 2D membrane can be cut into ID structures, so- 
called graphene nanoribbons (GNRs), which exhibit dif¬ 
ferent transport properties, depending on their termina¬ 
tion. Armchair terminated GNRs are usually gapped 
and therefore insulating. By virtue of their band g ap, 
they can be used to create field-effect transistors.®^ 
On the other hand, zigzag GNRs (ZGNRs) show local¬ 
ized edge states that may be spin-polarized.® Although 
armchair-type GNRs have been successfully synthesised 
using bottom up approaches,®® ZGNRs still remain elu¬ 
sive. Recently, patterned graphene with zigzag edge^^^ 
and GNRs with mixed armchair and zigzag terminations 
extending through a few lattice constant^ have been ex¬ 
perimentally realized. Despite the fact that the edges are 
not completely of zigzag type, they turned out to be of 
sufficient quality to confir m the prediction that the edge 
states become magnetized) 

Besides the geometrical confinement, another research 
area that has attracted much attention recently is the 
study of elastic deformations in graphene. Interest in 
this topic originated mainly from the theoretical predic¬ 
tion that strain may couple to the Dirac fermions as a 
pseudo-magnetic field (a magnetic field that preserves 
time-reversal symmetry). The subject was initially stud¬ 
ied in the light of deformations in carbon nanotubes.^ 
After the rise of graphene, this research direction grew 
in prominence by the vision of using strain as a way 
to tune graphene’s properties and use it in develop¬ 
ing an all-graphene electronics.^ The pursuit of strain 
engineering^^ was pioneered by the experimental obser¬ 
vation of “pseudo”-Landau levels in strained graphene,tI2l 
and has been recently corroborated by fascinating exam¬ 


ples of graphene spirals. 

In this paper, we study graphene systems that are 
both geometrically confined and strained, thus combin¬ 
ing the two research areas through a specific example: 
in-plane bent GNRs. These systems have been theo¬ 
retically investigated using a model based on density- 
functional theory.l^In addition, they have been proposed 
as a graphene geometry where strain couples as a uni¬ 
form pseudo-magnetic field.l^ Recently, these systems 
have been experimentally realized by pushing a GNR 
with the tip of a scanning tunneling microscope.!^ Al¬ 
though the experimentally synthesized bent samples are 
armchair terminated, here we concentrate on in-plane 
bent ZGNRs and study the dependence of the electronic 
behavior on the bending angle. We furthermore inves¬ 
tigate the dependence of the electronic structure on the 
type of bending. Our studies complement recent investi¬ 
gations of the mechanical properties of these systems.!^ 

As a main result, we find that the bending leads to 
an increased dispersion in the otherwise almost flat edge 
states. The bending breaks the symmetry between the 
inner and the outer edges, causing an effective compres¬ 
sion of the inside and elongation of the outside edge. To 
distinguish the contribution of each edge state to the dis¬ 
persion, we compare our findings to straight ZGNRs un¬ 
der compressive or tensile strain. Our results show that 
by tuning the bending angle, the edge states become dis¬ 
persive and hence delocalized. 

From a more abstract perspective, we can view the de¬ 
formations in graphene in terms of Japanese paper-art. 
Within this analogy, straight GNRs emerge as the art of 
paper cutting graphene. On the other hand, origami, the 
traditional Japanese art of paper folding, is connected to 
the study of strain in graphene. These two come together 
in graphene kirigami^ in which cutting and folding are 
combined. Here, the cutting refers to the specific ter¬ 
mination of the GNR, as well as to the fact that the 
hexagonal unit cells are empty (cut), and can thus be 
deformed in a variety of ways. A bent GNR is a very spe¬ 
cific and not so complicated type of graphene kirigami. 
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but precisely due to its relative simplicity, it is possible to 
study its electronic properties in depth. For this reason, 
the system is a good probe to understand how the elec¬ 
tronic behavior arising from confinement and termination 
is affected by strain. Therefore, we may generalize the 
notion of kirigami to more complicated graphene nanos¬ 
tructures, and apply a similar approach to understand 
their electronic properties. Knowing what is to be ex¬ 
pected in this simple case, may help us understand more 
complicated situations. 

This paper is organised as follows. In section |TI| we 
introduce a tight-binding model with exponentially de¬ 
caying hopping and overlap parameters, that we argue is 
suitable to study confined strained graphene systems. To 
the best of our knowledge, this particular tight-binding 
model has not been used previously to study GNRs, 
but turns out to capture all the relevant features of the 
band structure. We then introduce two types of bend¬ 
ing, which allow us to optimize the computational cost. 
In section |III[ we apply our model to study the effects 
of bending on the edge state and find that their localiza¬ 
tion can be tuned by bending. Conclusions are provided 
in section lYl 


II. A MINIMAL TIGHT-BINDING MODEL FOR 
BENT GNRS 

A. Three-parameter tight-binding model for 
strained confined graphene systems 

The electronic structure of graphene is usually derived 
using a tight-binding model with one Pz-OThital per site. 
If we assume a graphene system with n sites positioned 
at Vi, the single-electron wavefunction is given by 

n 

W) = . ( 1 ) 

i=l 

Here, are in general site-dependent basis states, 

which are assumed to be normalized. The vector c = 
(ci,... ,Cn)^ thus completely specifies the electron state. 
The Schrodinger equation can then be reduced to the 
n X n matrix equation 

{S£ + T)c = ESc, (2) 

where E is the energy associated with the state speci¬ 
fied by c. Here, we have split the Hamiltonian matrix 
H", the elements of which are given by Hij = {(j)i\H\(j)j)^ 
into the so-called hopping matrix (T), the diagonal on¬ 
site energy matrix (5), and the overlap matrix (S'), such 
that H = S£ + T. The elements of the overlap ma¬ 
trix are given by Sij = (0^10^). The matrix £ is diago¬ 
nal, with the elements corresponding to on-site energies, 
£ \(l)j) = ej |0j), as also defined in Ref. |25l Note that 
a more standard convention defines the on-site energy 
as the expectation value of the energy in a certain state 
and the hopping matrix as non-diagonal elements of the 


Hamiltonian matrix. However, the convention used here 
allows us to treat the on-site energy as a simple shift in 
E if the on-site energy is the same for each state. 

In general, we now have n(n + 1) parameters, the ele¬ 
ments of the matrices. In tight-binding, these parameters 
may be found by fitting to a reference calculation, rather 
than calculating them explicitly as integrals over basis 
functions. However, a model with n{n + 1) parameters 
is impossible to fit when n is not very small. Therefore, 
additional assumptions are made in order to reduce the 
parameter space. In graphene, translational symmetry 
allows one to use periodic boundary conditions. Since 
there is no longer a difference between individual sites, 
the on-site energy, hopping, and overlap parameters be¬ 
come site independent. 

A common procedure is to consider a two-parameter 
model that only takes nearest-neighbor (NN) and next- 
nearest-neighbor (NNN) hopping into account, and as¬ 
sumes orthogonal basis states. In this case, the site- 
independent on-site energy eo is left unspecified, as it 
leaves the eigenvectors invariant and produces only an 
absolute shift in the spectrum.lH However, we are inter¬ 
ested in a model for the graphene system that can de¬ 
scribe a bent GNR. For such a model, we have to specify 
the dependence of the hopping and overlap parameters on 
the distance, and, at the same time, the parameters of the 
model should not change when the system is geometri¬ 
cally confined, e.g. when graphene is confined to a GNR. 
This last condition would allow us to fit the parameters 
to a graphene reference calculation and not to a reference 
calculation for the specific GNR we study. We find that 
instead of the usual convention, a non-orthogonal model 
better satisfies these two conditions. First, we introduce 
the model and later argue why it compares positively to 
an orthogonal model. 

The tight-binding model we use is based on non- 
orthogonal site-independent basis states, which in real 
space are given by (r|0^) = 0(r — r^). Next to that, 
we assume that the hopping and overlap parameters be¬ 
tween these states are such that Uj = t{Tj — r^) and 
Sij = s{rj — Vi) where are exponentially decaying 
functions, given by 
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Here, a is the NN distance of graphene and to and sq 
are the values of the NN-hopping and overlap parameter, 
respectively. Note that the on-site hopping parameter 
is zero and that the overlap of an orbital with itself is 
one. The dimensionless constant k determines the fall-off 
rate of the hopping. Although this procedure introduces 
a discontinuity in the overlap that cannot be physically 
realistic, we will assume that the strain sizes are small 
enough, such that this effect can be neglected. We fur¬ 
ther assume that the hopping and overlap parameters 


3 


are proportional to each other, which implies that the 
parameter n is the same for both. 

This model satisfies the first condition we mentioned, 
a dependence of the hopping and overlap parameters on 
the distance, better than an orthogonal model. This can 
be seen by noting that in studies of strained graphene, 
exponentially decaying functions have been used fo r pa- 
rameters corresponding to orthogonal basis states.l^^^ 
However, efforts to reproduce the asymmetric band struc¬ 
ture of graphene using up to 20 fitted hoppings have re¬ 
sulted in subsequent parameters sometimes having op¬ 
posite signs and clearly not following a trend that can 
be described with an exponential decay.^^ On the other 
hand, if we relax the orthogonality condition, hopping 
and overlap are approximately exponentially decaying.^ 
When overlap is ignored in our parametrization (sq = 0), 
the model would be reduced to the one used in Ref. [26l 
Such a model does not reproduce the correct particle-hole 
asymmetry. Nevertheless, for low energies the overlap 
becomes less important and it would yield a good esti¬ 
mate of the spectrum. Orthogonal models which involve 
a non-exponential dependence on distance have also been 
used. Ref. [m for instance, introduces a separate linear 
dependence for both the NN and NNN hopping. One 
reason why this model is disadvantageous is that it has 
four fitting parameters instead of three, as in our case. 

An even more important reason for adopting the non- 
orthogonal approach is that these parameters are less 
dependent on the specific confinement than orthogonal 
parameters, thus better satisfying the second condition. 
To understand this, we note that in a quantum-confined 
graphene system we cannot expect all the hopping pa¬ 
rameters to have the same value as the bulk parameters, 
since now the edge needs to be taken into account. For or¬ 
thogonal states this is due, in part, to the fact that these 
states are a linear combination of P; 2 -orbitals obtained us¬ 
ing an orthogonalization scheme, like the Lowdin one.^^ 
These states are not the same on the edge and in the 
bulk, which also results in a difference of on-site energy 
and hopping between bulk and edge. Therefore, it is 
more realistic to assume non-orthogonal basis states for 
the tight-binding model. This allows us to get the param¬ 
eters from fitting to a graphene reference calculation and 
then apply it to the specific confined structure in which 
we are interested. A model based on nonorthogonal-basis 
states would be more universal than an orthogonal one 
for that reason. In Ref. mi an orthogonal tight-binding 
model is used and indeed we see that different hopping 
values are assumed for different GNRs: NNN hopping is 
zero for AGNRs and non-zero for ZGNRs. A more pre¬ 
cise way to treat the edge effect re quires the introduction 
of a different hopping at the edge.^^^^ However, for the 
sake of simplicity, we neglect this effect here. 

We have argued that the parameters of the model can 
be obtained by fitting to a reference graphene spectrum. 
In the periodic graphene case, Bloch’s theorem is used to 
reduce Eq. (§ to a 2 X 2 matrix equation, with wave func¬ 
tions labeled by the wavevector k in the Brillouin zone 



FIG. 1. (Color online.) Plots of the dispersion relation for 
graphene along the line connecting the T — M — K — T points 
of the Brillouin zone. The green-dashed curve corresponds 
to the two-parameter orthogonal dispersion of Ref. [U where 
the NN and NNN hopping parameters are t = —3.00236 eV 
and t' — 0.20509 eV, and eo has been chosen such that the 
K points are at zero energy. The blue-dashed curve depicts 
the orbital dispersion with exponentially decaying hopping 
and overlap parameter, described by Eq. Q, with values 
to = —2.8 eV, So = 0.2 , = 2.6, and eo = —1.28 eV, chosen 

such that the zero energy is at the K points. The red-dashed 
curve corresponds to an orthogonal-basis dispersion taking 
into account the first 15 hoppings,^ and the black-solid line 
corresponds to DFT calculations made using the Quantum- 
Wise software and a Hiickel type basis set.^ 


of graphene. In that case, the solution of this equation 
is equivalent to the one described in Ref. |29l By fit¬ 
ting to a reference first-principle spectrum, we find that 
to = —2.8 eV, So = 0.2, and n = 2.6 gives a reasonable 
match, which is also not very far off from the parame¬ 
ters used in Ref. [29l Although a more elaborated fitting 
method would allow us to find parameters that reproduce 
the reference spectrum more closely, we settle with these 
because we are mostly interested in global features and 
not in extremely precise quantitative results. 

The dispersion of graphene along a line connecting 
high-symmetry points of the Brillouin zone is shown in 
Fig-lU In this figure, different graphene dispersions ob¬ 
tained from different models are compared. One can ob¬ 
serve the results obtained from our three-parameter non- 
orthogonal model (blue-dashed line), the two-parameter 
orthogonal model of Ref. [1] (green-dashed line), and an 
orthogonal model where the first 15 hopping parameters 
of Ref. [23 are used (red-dashed line). The figure also de¬ 
picts the energy dispersion from a first-principle calcula¬ 
tion of graphene that was made using the QuantumWise 
software (black-solid hne).I^Erom the figure, we can ob¬ 
serve that the 15 parameter orthogonal basis model re¬ 
produces very well the dispersion relation obtained by 
first-principle calculations. The two-parameter orthogo¬ 
nal and three-parameter non-orthogonal models capture 
the essential features, but differ markedly at the M point 
for the chosen parameters. This is not surprising, as it 
has been shown that the behaviour around the M point 
is strongly influenced by higher-order hoppings 
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FIG. 2. (Color online.) Parameters characterizing the 
straight and the strained GNR. The parameters L and W 
characterize the straight ribbon, while the inner and outer 
radius Rin and Rout specify perfectly circular bending. Here, 
we also dehne the bending radius R along the center of the rib¬ 
bon, the width of the ribbon after bending W' = Rout — Rin, 
the length L' along the center of the bent ribbon, the length 
of a ID unit vector a along the center of the bent unit cell, 
the curvature Oid of a unit cell, the total curvature ^fuii of the 
ribbon and the bending parameter © = VF/2R. The region 
enclosed by the dotted line is the ID unit cell of the bent 
ribbon. 

B. Lattice-preserving bending 

To find a minimal model that can describe the ge¬ 
ometry of bent GNRs, we first introduce the concept of 
lattice-preserving bending. This type of deformation can 
be described by the parameters defined in Fig. We 
quantify the degree of bending using the dimensionless 
parameter 0, defined as 0 = hF/ 2 i?, where W is the 
width of the undistorted ribbon and R is the radius of the 
circular deformation. For W' « FF, this is approximately 
equal to the parameter used in Ref. [211 For straight rib¬ 
bons, we can define a ID unit cell with sites labelled by 
m, given by {r^}, and a ID lattice vector a. All sites 
can then be decomposed in + ia. for some num¬ 

ber £. This allows us to reduce the size of the matrices 
in Eq. ([2) using the ID Bloch’s theorem to 2N x 2N^ 
with N the number of dimer lines of the ribbon (num¬ 
ber of sublattice pairs in the unit cell, which is always 
even for ZGNRs), see for instance Ref. |T3| However, the 
ID translational symmetry that allows this procedure is 
broken after bending. A lattice-preserving bending is a 
type of bending deformation that still allows us to reduce 
the matrices in Eq. (§ to size 2N x 2N. This is possi¬ 
ble because a lattice-preserving bending F© satisfies the 
discrete rotational symmetry 

Fe(rj + a) = TZ-e^^Feiri), (4) 

where is the matrix that represents a clockwise 

rotation by angle Oid, and and a are the lattice sites 
and ID lattice vector of the straight ribbon, respectively. 


This symmetry can be seen as a type of modified pe¬ 
riodic boundary condition.!^ Because the Hamiltonian 
commutes with the rotation operator by an angle Oid, we 
write a ID Bloch-type wavefunction for a bent GNR in 
terms of a continuous quantum number k. In real space, 
Eq. 0 then assumes the form 

^e,fe(r) = ^ (5) 

£,m 

Here, k G [0, 27r], m runs over the atoms in the bent unit 
cell = F 0 (r^), and i runs over the number of unit 
cells in the ribbon. The vector = (cf..., 
therefore completely determines the electron state for a 
certain wavevector k and bending parameter 0. Namely, 
the components Cj of Eq. M are given by Cj = 
with j related to m such that rj = + ia. Erom 

the time-independent Schrodinger equation ([^, we can 
derive a matrix equation for the vector of orbital compo¬ 
nents C®’^, 

(5'©’^)-1t©’^c©’^ = (Ef - eo)c©’^ (6) 

Here, and Sk are 2N x 2N matrices with components 

QQ,k _ \ ^ iki c}('T? 

£ 

rpQ^k Dk£+('T> 

dmn y'<'-£eid^m 

£ 

where t and s are defined as in Eq. ® and D© is the 
spectrum of the eigenstates. In our calculations, we use 
the values for sq and k derived from graphene. The on¬ 
site energy is set to zero, giving a Eermi level close to, 
but not exactly at zero. After the calculation, the spectra 
are shifted by an amount eo to place the Eermi level at 
zero. As can be seen from Eq. 0 , the dispersion scales 
linearly with to when the scale is normalized around the 
Eermi level, and we can thus calculate the dispersion in 
terms of to without having to explicitly specify its value. 
The tight-binding model using non-orthogonal basis and 
exponentially decaying hopping and overlap in combina¬ 
tion with lattice-preserving bending may be used as a 
minimal model to study bent GNRs because it only re¬ 
quires three parameters and equations with matrices of 
size 2N X 2N. 


C. Two types of bending 

A realistic geometry for a bent GNR may be extracted 
from a molecular dynamics simulation, where bending 
affects both bond lengths and bond angles. The exact 
type of bending then depends on the ratio of the spring 
constants of the respective deformations. Erom previous 
work, it is known that the bond length in th e grap hene 
lattice is much stiffer than the bond angle.^^^^ This 
observation prompts us to explore a limiting scenario. 


— r 


nd\ 


^fld 
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FIG. 3. (Color online.) A bent ZGNR with N = 5 and 
© = 0.15 for bond length-preserving (red dots) and width 
preserving (blue dots). The width of the ribbon after bond 
length-preserving bending is Wb and after width-preserving 
bending is 

where bending is completely absorbed in bond-angle dis¬ 
tortions, and which we call hondlength-preserving bend¬ 
ing. In addition, we consider a distortion which we call 
width-preserving bending, where the atomic positions are 
rotated around a concentric point. The width-preserving 
bending is the same deformation as has been used in 
Ref. |22l Notice that the bond length-preserving bend¬ 
ing obeys the rules of graphene kirigami, since the paper 
can be folded (bond-angle deformations), but it cannot 
be strained (bond-length deformations). The fixing of 
the bond lengths in the bond length-preserving bending 
leaves the NN hopping unchanged, so that any perturba¬ 
tion in the electronic structure can mainly be ascribed to 
modifications of the NNN hopping. In contrast, bond 
lengths are allowed to change in the width-preserving 
bending scheme, so it may be expected that the changes 
in the dispersion are mainly due to changes in NN hop¬ 
pings. Comparing the effects of these two types of bend¬ 
ing on the spectrum, therefore, allows us to decouple the 
effects of NN and NNN distortions. 

Both bending deformations are depicted in Fig.[^ We 
can explicitly describe the width-preserving bending by 
the deformation function 

= ( 8 ) 

with r = (rx^ry). This deformation assumes that the 
ribbon is positioned such that the middle of the GNR is 
on the x-axis and the ribbon lies in the x^-plane. Hence, 
the y coordinate of the undeformed site is in the inter¬ 
val [—kF/2, kF/2]. One can easily verify that this bend¬ 
ing satisfies the definition of a lattice-preserving bending 
F^(r + ai(^, 0) = 7^_0^^F^(r, 0). This deformation is a 
width-preserving bending in the sense that the distances 
between sites in the direction along the width of the GNR 
remain unchanged. Another feature of this bending is 
that the strain in the direction along the ribbon width 
increases linearly from the inner to the outer edge. This, 
in conjunction with the fact that the bending considered 
here equally compresses on the inside as it stretches out¬ 
side, yields a line of zero stress exactly in the middle of 
the ribbon. 


It is not straightforward to give a closed formula for 
the bond length-preserving bending. However, we can 
construct the profile of the deformation by applying 
F^^(Ri,0) on specific ribbon sites recursively, see 
Appendix. The bondlength-preserving bending is sim¬ 
ilar to the width-preserving one, but has a non-linear 
strain profile from the bottom to the top of the ribbon. 
At the inner edge, the ribbon experiences not only lon¬ 
gitudinal compressive strain, but also transverse tensile 
strain. At the outer edge, on the other hand, a com¬ 
pressive transverse strain is present. It is also important 
to note that the total width becomes reduced, as can 
be seen in Fig. This reduction of width needs to be 
taken into account when comparing effects of the bond 
length-preserving with the width-preserving bending. As 
a consequence of the reduction of width, the longitudinal 
strains at the inner (e^^) and outer edge (e^ut) are not 
identical for the two types of bending. 


III. RESULTS: TUNABLE EDGE STATE 
DISPERSION 

We have calculated the dispersion relation for bent 
ZGNRs by solving Eq. ^ numerically both for width¬ 
preserving and for bond length-preserving bending. In 
Fig. i the dispersion relation for different values of 
the bending parameter is depicted. Since we argued 
that bending introduces a profile of elastic deformation 
with effective compressive strain on the inside and ten¬ 
sile strain on the outside, it is useful to compare it to 
the effects of a uniform longitudinal strain e, defined as 
e = AL/L, where AL is the length deformation intro¬ 
duced by the strain, and L is the length of the unde¬ 
formed nanoribbon. Fig. depicts the effect of positive 
(tensile) and negative (compressive) longitudinal strain 
on a A = 4 ZGNR subjected to a width-preserving uni¬ 
form strain deformation. We can see that the energy of 
the edge states increases (decreases) for negative (pos¬ 
itive) strain. When we compare the two cases with a 
ribbon bent using width-preserving bending, we observe 
that the energy increase in the edge state that experi¬ 
ences compression is roughly equal to the energy increase 
in both edge states of a longitudinally compressed rib¬ 
bon. Similarly, we find a good agreement for the outer 
edge state with both edge states of a ribbon experiencing 
tensile strain. These observations indicate that the dis¬ 
persion of ribbons bent by 0 is quantitatively related to 
the dispersion of a uniformly strained ribbon with strain 
e = ±0, a result consistent with Ref. nn 

Plotting the wavefunctions of the edge states confirms 
that the low-energy state resides on the outside, as shown 
in Fig.[^ Here, the orbital components of the eigenstates 
of the edge states, Cj, as defined in Eq. (§, are plotted 
for increasing 0. The width-preserving bending scheme 
was used in generating the plots. Eirst, we note that the 
edge states are localized on one sublattice at both edges, 
forming a symmetric and antisymmetric combination of 
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states localized on either edge. The states are almost 
degenerate, which would allow us to form ortho normal 
combinations that are still eigenstates of the Hamilto¬ 
nian with the same eigenenergy. In fact, since interac¬ 
tion effects arising from the Coulomb repulsion are not 
accounted for, we may expect these effects to favor a dif¬ 
ferent combination in the two-dimensional Hilbert space 
of eigenstates. Intuitively, the effect of electron-electron 
repulsion should be to split the symmetric and antisym¬ 
metric states into two states that are localized on their 
respective edges, both singly occupied. 

With increasing bending, we observe that the nearly 
degenerate states that initially reside on both edges in our 
model transform into a high-energy state localized on the 
inner edge and a low-energy state localized on the outer 
edge. It is interesting to note that this already occurs 
for the very small bending parameter of 0 = 10“^, indi¬ 
cating that for this strength of bending, the symmetric 



FIG. 4. (Color online.) Dispersion relation and DOS for a 
N = 4 ZGNR as a function of the bending parameter © for 
width-preserving bending (blue) and bond length-preserving 
bending (red). The thin-black line corresponds to a straight 
ribbon. On the left side, we show the spectrum over the 
complete BZ, for 0 varying from 0 at the top to 0.15 at 
the bottom panel, with steps of 0.05. On the right side, we 
zoom in on the edge state with k ranging from 27r/3 to 47r/3 
and depict the lattice of the GNR. All plots have the same 
scale as shown in the bottom. Calculations were made using 
non-orthogonal parameters with exponential decay, given by 
Eq. with So = 0.2 , = 2.6, and eo such that the Fermi 

energy (dotted line) of the straight ribbon lies at zero. The 
DOS is calculated using a Lorentzian broadening with a width 
of 0.03 eV (DOS in arbitrary units). 


and antisymmetric states mix in order to form the states 
localized on a single edge, energetically more favorable. 
We also find a significant dependence of the localization 
length of both edge states on the momentum k. When 
we plot, for example, the edge states for a wave vector of 
k = Ttt/S, the wave function appears to spread more into 
the bulk of the ribbon than for the value /c = tt, as shown 
in Fig.|^ Although not shown here, the results for bond 
length-preserving bending show that the edge state for 
k = Ttt/S is also less localized than for k = it. However, 
for the same degree of bending, the effect is much less 
pronounced than for width-preserving bending. 

Another striking observation is that the two edge states 
do not only split but also develop opposite curvature, as 
shown in Fig. The top band is curved upward, but 
at its center a small downward curvature develops, such 
that there is a local maximum at /c = tt, whereas the 
opposite occurs for the lower band. This is in contrast 
with what we observe for positive or negative uniform 
strain in Fig. In that case, the edge states are only 
shifted, but retain the same dispersion as in the strain- 
free ribbon. 

A minimal model that captures this behaviour, and in 
particular fits the dispersion of the edge states around 
the point /c = tt, is a tight-binding model of a ID chain 
of sites with a single NN hopping and an on-site 
energy Here, the superscripts refer to the higher- 

energy band and lower-energy band, which are localized 
on the inner and outer edge, respectively. The effective 
dispersion obtained from the ID NN tight-binding model 
reads 

cos(fc). (9) 

Inspection of the zoomed in panels of Fig. [^suggest that 
this effective model can describe the shape of the bands 
in the region around k = tt reasonably well. A positive 
or negative relates to the dispersion that exhibits, 
respectively, an upwards or downwards curvature around 
momentum k = tt. 

Before we compare this effective model quantitatively 
with the tight-binding results, we need to mention the 
effect of the width of the ribbon on the edge states. As 
a ribbon becomes narrower, the edge state localized on 
one side with k closer to k = tt starts to hybridize with 
the edge state localized on the other edge. On the other 
hand, when one starts bending a ribbon the edge states 
start moving closer in energy to the bulk states. This 
can be seen in Fig. After a certain bending, the va¬ 
lence band maximum hybridizes with the lower-energy 
edge states, as well as the conduction band minimum hy¬ 
bridizes with the higher-energy edge states. Since wider 
ribbons have a smaller bulk band gap, these effects are 
more pronounced. These effects are shown in Fig. 
where we plot the lower-energy edge state for three dif¬ 
ferent widths of the ribbon, = 4, 14, and 30, using 
the same bending parameter, 0 = 0.1, and /c—value 
k = Gtt/S. We observe that the two edge states of the 
N = 4 ribbon hybridize with each other, and are there- 
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e=-0.1 


e=0.1 



e=-0.1, £=0.1, 0=0.1 



6 2it/3 4it/3 2tt 


FIG. 5. (Color online.) Dispersion relation for a iV = 4 ZGNR for uniform strain e = —0.1 (left panel, red line) and e = 0.1 
(middle panel, blue line). The thinner-black line in the left and middle picture corresponds to the straight ribbon. The right 
panel shows the dispersions for uniform strain e = —0.1 (red), e = 0.1 (blue), and width-preserving bending with © = 0.1 
(black). Calculations were performed using non-orthogonal parameters with exponential decay, given by Eq. § with So = 0.2 
, Hi — 2.6, and eo such that the Fermi energy (dotted line) of the straight ribbon lies at zero. 


fore not localized anymore. The edge states of the ribbon 
with = 30 also hybridize, but instead with bulk states, 
and are not localized anymore either. The ribbon with 
A/" = 14, however, still shows localized edge states for the 
same regime of parameters. These two opposite effects 
make the comparison between different ribbon sizes very 
intricate. We have chosen to analyze the N = lA ribbon 
in more detail because this one has the optimal width 
to avoid spurious hybridization effects of the first or sec¬ 
ond kind. Our observations are expected to hold also for 
ribbons of different width, if care is taken to account for 
these hybridization effects. 

We fit the parameters of the effective ID dispersion 
of Eq. to the tight-binding calculations for a ribbon of 
width N = 14. In Fig. we plot the fitted parameters 
for different values of the bending 0. We observe that 
both the lower- and the higher-energy edge states start 
out with the same positive hopping parameter. Interest¬ 
ingly, in both bending schemes, crosses zero, implying 
that for a certain bending parameter the band becomes 
dispersionless. This is an important observation because 
many-body effects can be expected to become even more 
relevant for that bending parameter. 

By comparing how the parameters change with respect 
to the type of bending used we can identify whether the 
NN or the NNN hopping is more important. The effective 
parameters for the state on the inner edge decrease for 
both types of bending. However, for the outer edge the 
effective parameters increase for bondlength-preserving 
bending, but decrease for width-preserving bending. The 
main difference between the two bending methods is that 
in the width-preserving bending also the NN distance is 
modified. Therefore, we can conclude that for the outer- 
edge state NN effects are more important than for the 
inner-edge. General behavior of the inner-edge state, 
however, can be captured by only considering the effect 
of the NNN hopping. If we compare the effective param¬ 
eters for the inner-edge state between the two bending 
methods in more detail, we observe that the effective pa¬ 


rameters for bond length-preserving bending show a lin¬ 
ear dependence on 0, while this dependence for width¬ 
preserving bending is more complicated. One reason for 
this behavior could be the fact that the width- and bond 
length-preserving bending produce a small difference in 
strain on the edges (e^^, Cout)- To check whether this 
can account for the difference, we also plot the effective 
parameters as a function of the strain (smaller plots in 
Fig.[^. We can clearly see that the general behavior does 
not change. Therefore, the difference should be sought in 
effects of the NN hopping. Changes in the NN distance 
influences the hybridization between the opposite edges 
and the hybridization of the edge state with bulk states. 
These effects might explain why the effective parameters 
of width-preserving bending exhibit a nonlinear depen¬ 
dence on the bending. Furthermore, the effect of the per¬ 
turbation of the NN distance also depends on the width of 
the ribbon, which additionally complicates the problem. 
Because of all this, in the following we focus only on the 
effective parameters of bond-length-preserving bending. 

For bond length-preserving bending (see plots in red in 
Fig.[^, the effective hopping at the inner edge (higher-F^) 
linearly decreases and changes sign, whereas the hopping 
at the outer edge (lower-F^) linearly increases. We could 
try to understand this behaviour by assuming a perfectly 
localized edge state. The inner edge experiences a neg¬ 
ative strain, so the hopping becomes more negative and 
the ID dispersion would curve downwards. This indeed 
corresponds to what we observe in Fig. On the same 
token, the hopping at the outer edge should decrease, be¬ 
cause the distances between the lattice sites increase, and 
therefore a flat band should develop. However, the oppo¬ 
site behaviour is visible in Fig.[^ This can be understood 
by noting that the changes due to bending at the outer 
edge are determined by the weight of the wavefunction 
on sites closer to the bulk. This is because these sites are 
closer to each other and therefore contribute more to the 
energy. This together with the fact that sites close to the 
bulk have a sizeable weight implies that our assumption 








































of the localized states does not apply. The fact that the 
edge state becomes less localized as the momentum moves 
further away from k = ir is crucial here. This enhances 
the effect that can already be seen for straight ribbons, 
where the edge states are dispersive at the momenta away 
from /c = TT, and causes an increasing positive effective 
hopping. 

In conclusion, we can understand the behaviour as a 
competition between two effects due to NNN hopping 
and strain: 

1. An effective positive hopping for increasing neg¬ 
ative strain because of the increasing delocalized 
nature of the edge state as the momentum moves 
further away from k = n. 

2. An effective negative hopping for increasing nega¬ 
tive strain because the edge state is localized. 


Color code for complex phase cp 



FIG. 6. (Color online.) Edge states in the real space. Compo¬ 
nents Q,rn of edge states for a section of the ribbon at /c = tt 
and k = Ttt/S are mapped to the corresponding points of 
the bent GNR for different values of bending-parameter ©. 
The width-preserving bending scheme was used here in com¬ 
bination with our standard orbital hopping parameters, q m 
are related to the eigenvector through the definition Eq. ra 
and satisfy Eq. The coefficients Q,rn are complex num¬ 
bers that are depicted in two ways. The first is by dots of 
which the diameter is proportional to the absolute value of 
Q,rn and the color corresponds to the phase, as indicated by 
the color code. Additionally the coefficients Q,rn are repre¬ 
sented by a vector in the complex plane (see small black lines 
at the center of the dots). The phase is chosen such that the 
lattice site at the left bottom of the picture has phase zero. 
The top (bottom) rows of the /c = tt and k = Ttt/S panels 
correspond to the high-energy (low-energy) edge states. The 
ribbon has width iV = 6. 


For the edge state localized on the outer edge, the first 
effect is always dominant and becomes even more rele¬ 
vant after bending. For the inner edge, the second effect 
overcomes the first after a certain bending parameter. 
This is the reason why the dispersion of the inner edge 
has to go through a point at which it is dispersionless. 
This also clarifies our earlier observation that the outer 
edge state is more sensitive to changes in the NN hop¬ 
ping. The outer edges are more delocalized, and therefore 
more sensitive to the effects of the NN hoppings. 



FIG. 7. (Golor online.) Real space depiction of lower energy 
edge states for ribbons of with A = 4, A = 14, A = 30. 
Gonstruction is the same as explained in Fig.[^ For all three 
ribbons the bending parameter is 0 = 0.1 and k = Gtt/S. The 
A = 4 ribbon shows no clear edge states because the states are 
hybridized across the entire ribbon, while the A = 30 ribbon 
shows localized edge states, which are, however, hybridized 
with ones in the bulk . 
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FIG. 8. (Color online.) Best fit values of the effective ID chain parameters Cq (upper rectangular panel) and Cq (lower 
rectangular panel) versus the bending parameter 0. These parameters are dehned in Eq. The htting has been performed 
for width-preserving bending (blue dots) and bond length-preserving bending (red dots). The right picture in each panel shows 
these parameters with respect to the strain on the inner edge Cin for the higher-energy edge state and with respect to the 
strain on the outer edge Cout for the lower-energy edge states. The relation between e^n, Cont and 0 is explained in the text. 
The htting was done for a ZGNR of width N = 14, and is based on data points chosen in the region around /c = tt given by 
k E [2.41, 3.86]. Error bars are obtained from the standard deviation between the htted spectra and the numerics on the lattice. 
All calculations are performed using our standard set of exponentially decaying orbital hopping parameters. 


IV. CONCLUSIONS 

We show here that a tight-binding model with ex¬ 
ponentially decaying hopping and overlap can be used 
as a minimal model with three parameters to study a 
graphene-based system that is both geometrically con- 
hned and strained. To obtain geometries of bent nanorib¬ 
bons that serve as the input of the tight-binding model, 
we develop two types of bending, bond length-preserving 
and width-preserving. We would like to point out that 
bond-length preserving bending geometry, generated us¬ 
ing a recursive algorithm, shows a particularly strong 
analogy with the Japanese art of kirigami. Both types 
of bending are lattice-preserving, causing the resulting 
geometry to exhibit rotational symmetry (the unit cell 
is rotated by Oid to generate the entire bent GNR), and 
therefore allowing us to reduce the tight-binding model to 
the numerically inexpensive problem of solving a matrix 
equation with 2N x 2N matrices, with 2N the number 
of sites in the unit cell of the GNR. The different types 
of bending allow us to decouple the effects of perturba¬ 
tions of the NN and NNN parameters of the tight-binding 
model. 

We have investigated the qualitative features of the 
dispersion relation upon bending. Our calculations show 
that bending leads to nontrivial effects on the edge states 
of ZGNRs, resulting from the broken symmetry between 


the top and bottom edges. We observe that both width¬ 
preserving and bond length-preserving bending predict a 
splitting of the two edge states (without considering in¬ 
teractions). A lower-energy edge state localizes on the 
outer edge and a higher-energy edge state on the inner 
edge. In fact, there is an emergent band structure around 
the point k = tt of the edge states that can be fitted to 
the tight-binding dispersion of a ID chain with an effec¬ 
tive hopping and on-site energy parameter. The higher- 
energy edge state has an effective hopping parameter that 
changes sign as the bending is cranked up from 0 = 0.11 
to 0 = 0.17, with the exact value where the effective hop¬ 
ping vanishes depending on the type of bending. Hence, 
there is a critical degree of bending at which the band 
is effectively fiat and interaction effects are expected to 
become prominent. Since the charge carrier velocity is 
proportional to the slope of the dispersion, the degree of 
localization of the edge states can be tuned with bend¬ 
ing. By comparing the two bending methods, we can 
conclude that effects on the dispersion of the inner-edge 
state are dominated by changes in NNN hopping. For 
the outer edge state, changes in NN hopping also become 
important. The effects due to NN hopping changes, how¬ 
ever, are less universal and depend on width and bending 
method. The effects of the NNN hopping on the emergent 
band structure at the edges can be explained by a com¬ 
petition between the decreasing localization of the elec¬ 
tronic states with the momenta away from k = n and the 
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localized character of the edge state. A next step would 
be to include interaction effects, as these are important 
for edge states, especially when the bending gives rise to 
the flat bands. Furthermore, motivated by our work, it 
would be important to understand how bending would 
affect the magnetic polarization of the edge states de¬ 
tected recently.^ We hope that our results will stimulate 
further research in these directions. 
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APPENDIX 

Recursion formula for bond length-preserving 
bending of ZGNR 

We construct the bond length-preserving bending, 
for a ZGNR. First, we construct the bent ID unit cell. 
The orientation is chosen such that the first site in the 
bent unit cell is positioned at = (0, R — lF'/2). Note 
that we do not know W' and R yet, but they will be 
obtained using a recursive procedure outlined below. We 
can now recursively generate the next atoms in the de¬ 
formed ID unit cell using the following rule: 

/Id ^ f 6'id)^-(?id/2rm-i if * is even 

1 (kmli I + a)fm-i if * is odd, 


= |r;ili|cos(6»id/2) 
+ \/(a)^ - |rm-iPsin2(6>id/2). 


Here, is the unit vector in the direction of 

We still assume that the distance along the middle 
of the GNR remains unchanged, and therefore Oid = 
Qa'/(W/2). If we follow this recursion until r^^^, where 
N is the number of A sites in the ID unit cell, we have 
generated the deformed ID unit cell However, we 

started with defined in terms of the bent GNR width 
VF', which was unknown. We can now use the identity 
W' = |r 2 ^^| — |r'/^|, which is an equation with W' on 
both sides, to write out the recursion explicitly. How¬ 
ever, this is a rather involved equation. We can, on the 
other hand, easily find a good approximation iteratively 
for W'. We start with the assumption that W' « W. 
Then, after running the recursion, we calculate the W' 
of that ribbon. If it differs by more than a set test value 
from the previous recursion, we use that value of W' to 
generate a new unit cell. This iterative procedure runs 
until the test condition, that gives the minimal differ¬ 
ence between a new and old width, is satisfied. Note also 
that this deformation does not work for every 0, as for 
large enough bending the square root in the definition 
will become complex. This is understandable, as there 
should be a maximum bending at which the lattice sites 
on the outer edge of the ribbon are all separated by a. 
Once the bent unit cell is generated, the complete bent 
GNR is obtained by copying the unit cell through multi¬ 
ples of rotations by Oid. Thus, we can describe the bond 
length-preserving bending as 

Here, £ runs over the number of unit cells in the ribbon. 
We explicitly use that the lattice sites of a straight GNR 
can be described by a site in the ID unit cell plus a 
multiple of a, the lattice vector of the straight ribbon. 
One can show, using simple trigonometry, that each site 
now has 3 neighbors that are at a distance equal to a, as 
shown in Fig. Due to the construction, it is obvious 
that the rotational symmetry is satisfied and thus this is 
a lattice-preserving bending. 
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